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We present a determination of nucleon-nucleon scattering phase shifts for £ > 0. The S, P, D and F 
phase shifts for both the spin-triplet and spin-singlet channels are computed with lattice Quantum 
ChromoDynamics. For ^ > 0, this is the first lattice QCD calculation using the Lirscher finite- 
volume formalism. This required the design and implementation of novel lattice methods involving 
displaced sources and momentum-space cubic sinks. To demonstrate the utility of our approach, 
the calculations were performed in the SU(3)-flavor limit where the light quark masses have been 
tuned to the physical strange quark mass, corresponding to m.n-=mx~800 MeV. In this work, we 
have assumed that only the lowest partial waves contribute to each channel, ignoring the unphysical 
partial wave mixing that arises within the finite-volume formalism. This assumption is only valid for 
sufficiently low energies; we present evidence that it holds for our study using two different channels. 
Two spatial volumes of E ~ (3.5fm)^ and V ~ (4.6fm)^ were used. The finite-volume spectrum is 
extracted from the exponential falloff of the correlation functions. Said spectrum is mapped onto 
the infinite volume phase shifts using the generalization of the Lirscher formalism for two-nucleon 
systems. 

PACS numbers: ll.80.Et, 12.38.Gc, 13.75.Cs, 13.85.Dz 


I. INTRODUCTION 

Understanding low-energy nuclear physics directly from the underlying theory of strong interactions, Quantum 
ChromoDynamics (QCD), remains a primary goal of nuclear physicists. The motivation can be broadly separated 
into two categories: obtaining a quantitative description of basic nuclear physics directly from QCD and probing the 
limits of the Standard Model and its fundamental symmetries through precision low-energy experiments in nuclear 
environments. In both cases, there are substantial international experimental efforts planned or underway which 
require a quantitative understanding of QCD. 

The basic interactions of two nucleons (NN) and nuclei are well measured and have led to a variety of precise 
theoretical descriptions ranging from phenomenological models to effective field theories (EFT). These over-constrained 
NN interactions are then used to predict properties of light nuclei using a variety of methods such as multi-nucleon 
EFT [IH3] I harmonic oscillator based effective theory |3J , no-core shell model [HI [7] and Green’s Function Monte 

Carlo [an]. For light nuclei, the NN interactions dominate the nuclear structure, but the three-body nuclear force is 
necessary for accurate comparisons with the measured values [Hill]. A recent, exciting development is the use of 
lattice field theory to regularize the two- and three-nucleon EFT and predict properties of light nuclei, such as the 
recent computation of the Hoyle State [T7] . 

All of these impressive theoretical applications rely upon experimental input of the nuclear interactions. While 
this is achievable precisely for the low-energy NN interactions, there are very few constraints on the three and higher 
nucleon forces. Determining these interactions directly from QCD is a multifaceted problem. At the core of this 
challenge is the non-perturbative nature of QCD which requires a numerical approach at low energy. Lattice QCD 
(LQCD) is the discretized version of QCD in a finite Euclidean volume. It is the only known tool to compute QCD 
correlation functions in the infrared for which no uncontrolled approximations are necessary. Using LQCD combined 
with EFT (see for example Refs. [nmi]), we will be able to determine the elusive few-body nucleon interactions 
directly from QCD, relevant for example for the upcoming experiments at FRIB m, designed to study neutron-rich 
nuclei. We will also be able to compute the hyperon-nucleon interactions, which are extremely challenging to measure 
experimentally due to the rapid weak decay of the hyperons |16H19j . These calculations will be relevant for the 
experiments planned at the FAIR, JLab and J-PARC facilities. There are even studies of hyper-nuclei using heavy 
ion collisions at RHIC |20] and the LHC [TT] . 

LQCD is also necessary to compute one, two and few-body nuclear matrix elements such as the scalar matrix ele¬ 
ments needed for direct dark matter detection and the electroweak matrix elements which govern nuclear interactions 
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and decays. Two recent examples of these include the N —)■ Ntt parity violating process [22] and the parity conserving 
np —>■ dj rate [23] . It is known that, due to technical complications, calculations of matrix elements involving multi¬ 
particle states often receive Oil) corrections from the finite volume P^26j . In order to control these corrections, the 
two-particle phase shifts and their derivatives must be determined. See Ref. m for a very nice demonstration of this 
technology for the case of 717 —>■ tttt. 

Tremendous progress has been made in performing LQCD calculations of two-meson interactions [28l - l36| . These 
calculations use the Liischer formalism |37H45j to relate energy levels in a finite periodic volume to the infinite volume 
scattering phase shifts. More recently, this has been extended to include coupled channels such as the TrK-rjK 
system gillTI and the tttt-KK 1=1 channel jUj. In contrast, the NN system is much more challenging to study 
for a variety of reasons, see for example [4^ [50] . These calculations have been limited to ^-wave interactions and 
bound states [siHsni.^ The methodology for determining three-body forces from LQCD calculations is still being 
developed |55H55] . However, it is evident that to reliably extract three and four-body forces a precise determination 
of two-body scattering parameters is needed, including but not limited to £ = 0 partial waves. 

In this work, we present a calculation of higher partial wave scattering in the NN system. In particular, we have 
computed S,P,D and F partial waves in both isosinglet and isotriplet channels using the NN generalization [^ of 
Liischer’s formalism. Given the complex nature of this problem, this exploratory calculation was performed at the 
S'C/(3) flavor symmetric point with ~ 800 MeV, enabling us to explore the implementation of our new method 
and demonstrate its feasibility with relatively little investment of computing time. We have also simplified the 
Liischer formalism by ignoring mixing from higher partial waves contributing to a given cubic irrep. This mixing 
is kinematically suppressed at sufficiently low scattering energies, however, for the range of energies we explore the 
assumption that they do not contribute significantly may require further investigation. Some evidence that the 
assumption is valid can be obtained by comparing different cubic irreps coupling to the same partial wave, and is 
presented in Section III This work is an extension of a previous determination of the NN S'-wave interactions [53] on 


the same LQCD gauge configurations. 


II. IMPROVED TWO-NUCLEON INTERPOLATING FIELDS 


We have performed the calculations with the isotropic Wilson lattices generated by the JLab/W&M group, with a 
lattice spacing of 6 = 0.145(2) fm and two spatial extents of L/b = 24 and L/b = 32 (for more details see Ref. [57] 1. 
Since the volumes are cubic and the two-nucleon systems considered have zero total momentum, the spectra obtained 
are those of the irreducible representations (irreps) of the octahedral symmetry group Oh- 

The sink operators are defined as products of single nucleon operators in momentum space. The single nucleon 
operators were designed to have good overlap with the single nucleon ground states [701172] . Let R be an element of 
Oh- Let (k) be a nucleon operator with spin and isospin z-components rris-^ and mi^ and momentum k at time 
t. Due to the periodic boundary conditions, the free momenta satisfy k = 27rn/L, where n is an integer triplet. With 
these we can construct two nucleon operators with angular momentum Jmj and isospin Imj 


O 


Jmj 

Imi\S£. 


|k|)= E 


ms ,'m^ 
TTLs-^ ,ms2 
mi-^ 


Jmj ^Sms ^Imj 

R£Oh 


^ r,™,(Rk) vri((i?k)vrii(-i?k). 


( 1 ) 


where i and S denote the total orbital angular momentum and spin of the system and i?k is the unit vector in 
the Rk direction and the Y^-meiRk) are the standard spherical harmonic functions. The standard Clebsch-Gordan 
coefficients, C;^Z,’sms = {Jmj\£mi, Sms) project the operators to total spin S, angular momentum J, and isospin /. 
The infinite volume quantum numbers are not good quantum numbers in a finite volume. Operators with different 
angular momentum labels will mix due to the non-spherically symmetric finite spatial boundary conditions. To project 
the operators above to an operator that is in a row p of the irrep A of Oh^ we use the subduction coefficients, 
found in Ref. [73| 


mj 


( 2 ) 


^ There has been one exploratory study of higher partial wave NN interactions in Ref. m- This calculation used the so-called potential 
method m which suffers from additional systematics that have not yet been demonstrated to be under control. To date, there has 
been one quenched 7 = 2 tttt comparison m- 
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FIG. 1: Effective masses for the energy splitting, AE„ = — 2mjv, in lattice units for the second excited state in 

the spin triplet channel at L/b = 32, showing operators having different [JIS] labels (Eq. Q): J = 1,1 = 2, S = 1 (black), 
J = ?>,l = 2, S = 1 (blue), J = 1,1 = 0, S' = 1 (red), J = 1,1 = 0,S = 1, r = 0 (green). The dashed horizontal lines represent 
the energy levels of the nearest non-interacting two-nucleon states. 


While of course it is possible to bypass the construction of operators with [JIS] labels in Eq. Q and go directly to 
a set of operators belonging to definite cubic irreps, it is convenient to keep separate operators of the same cubic irrep 
having different [JiS] labels, as these operators in some cases have very different overlap with the various excited 
states of a particular cubic irrep. For example, in the spin singlet, channel, the second non-zero momentum shell 
of the non-interacting system contains two degenerate states, whose energies split once interactions are turned on. 
We find that the operators having £ = 0 labels exhibit good overlap with the lower of these two states, while 
operators with £ = 2 labels overlap well onto the higher state (see Fig. [^. 

Ideally, we would construct NN operators in momentum space at both the source and sink locations, as is done in 
two-meson calculations [Ml ES] • However, the computational cost of performing these calculations for two nucleons is 
orders of magnitude greater than for two mesons. For this reason, NN calculations typically are performed with local 
or volume sources for the nucleons at the source and then projected to definite momentum at the sink. 

The significant improvement made over previous works is the use of spatially displaced nucleon operators at the 
source, fV'^(£o,xo-l-r/2)iV'^(to,xo —r/2). By displacing the two nucleons at the source, we find a significant increase in 
the overlap of the operators onto the NN H)*" and irreps (~ S'-wave) as compared to the local operators (see Fig. 
left). Further, without such a displacement, for zero total momentum, the overlap of the local operators (r = 0) with 
the cubic irreps that contain the P, D and F waves are zero, prohibiting a determination of the spectrum in these 
irreps. 

Displaced operators give us further freedom in designing our sources to have good overlap with the desired states 
by choosing from various geometries for the displacements. For example, after fixing one nucleus to a single lattice 
point, (0,0,0), one may then calculate the set of pairs where the second nucleon is displaced in all possible ways 
along a single axis, r = |r|(0,0,1) (plus all cubic rotations), or along multiple axes, such as r = |r|(0,1,1) (plus 
all cubic rotations) and r = |r|(l,l,l) (plus all cubic rotations). We have named these geometries “face”, “edge”, 
and “corner”, respectively, and disregard more complicated geometries. Each collection of geometries may then be 
projected onto the desired cubic irrep as described above for the sink operators (Eq. Q), with the momentum vectors 
k replaced by the set of displacements vectors, r. Note, however, that this is only a partial projection because we 
use a reduced set of displacement vectors compared to the full lattice volume. Example effective masses for the three 
types of geometries are shown in Fig. (right). While the “face” sources have a reduced computational cost compared 
to “corner” and “edge” sources (7 inversions versus 9 (“corner”) and 13 (“edge”)), they have zero overlap with several 
channels of interest due to their simple geometrical structure. We have chosen to focus on “corner” sources for the 
remainder of this work, to balance good overlap with a large number of states with moderate computational cost. 

In Fig. we show effective mass plots for several operators in two different cubic irreps along with the resulting 
determination of the energy levels. We construct all sink operators that have free momenta |k|L/27r < VS. Note that 
the TJ channel has overlap with the physical ^P 2 , and channels, and the channel has overlap with the 
physical ^D 2 channel. Both channels have additional overlap with kinematically suppressed higher partial waves. The 
bands denote the combined statistical and systematic uncertainties obtained by performing one and two exponential 
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FIG. 2: Effective masses for the energy splitting, AE„ = 2\J+ Qn ~ in lattice units for the first excited state in the 

spin singlet channel at Ljb — 24. On the left, we show effective masses corresponding to different displacement lengths at 
the source: |r| = 0 (green), |r| = 1 (red), |r| = 3 (black), |r| = 5 (blue). On the right, we show effective masses corresponding 
to different geometric displacements at the source (see text): face (blue), corner (black), edge (red) (“face” sources calculated 
on a significantly smaller sample of configurations). The dashed lines represent the energy level of the nearest non-interacting 
two-nucleon state. 
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FIG. 3: Effective masses for the energy splittings, AEn = — 2mjv, in lattice units for spin triplet Tj and spin 

singlet T 2 channels. Several energy levels are shown (corresponding to the non-interacting shell labeled by n), and red circles 
(black squares) points correspond to L/b — 24 (32). Horizontal bands represent fits to the data (see text). Dashed (dot-dashed) 
lines represent the energy levels of the non-interacting two-nucleon systems. 


fits to the correlation functions over a range of time windows. The fit ranges displayed are representative of the times 
considered in the fits. For reference, the dashed lines correspond to the energy levels of a non-interacting two-nucleon 
system. The calculations were performed with 20 sources on 3822 configurations and 70 sources on 1018 configurations 
for the L/h = 24 and L/b = 32 volumes respectively. 


III. NN SCATTERING PHASE SHIFTS 

In general, due to the reduction of rotational symmetry in a finite volume, there is not a one-to-one correspondence 
between the finite volume spectrum and the infinite volume scattering amplitudes. For sufficiently low energies, one 
can expect higher partial waves to be kinematically suppressed. Ignoring partial wave mixing and .^ > 4 waves, the 
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FIG. 4: Shown are examples of the phase shift determination in several partial wave channels, as well as representative ERE 
fits. The dashed vertical lines indicate the momentum at which the t-channel cut occurs {q — m,r/2). In most panels we plot 
q^^~^^cot5e which is used to determine the parameters in the ERE. In the ^So channel, the faded points were not included 
in the ERE fit, but are displayed here to show consistency. In the upper right panel we also show the phase shift as a 
function of the lattice momenta. 


spectra of a number of cubic irreps satisfy the quantization condition 


gcot^A(g) 

47r 


_ / 2n , C4o(g'^) , C6o(g^) 

— COo(9 ) + CI!4,A-4-1" Q!6,A- q - 


( 3 ) 


where q is the on-shell relative momentum of the system, 5f^{q) is the scattering phase shift of the 

partial wave that primarily couples to the A irrep, o^^a are constants, Ref. m, and the are kinematic, non-linear 
functions that depend solely on the momentum and the volume 


{q ) 


L3 I L / 2^ (r2 - q2) • 

\ / rGZ3 ^ ^ ' 


( 4 ) 


Employing Eqs. ^ and ( [4|), we obtain the scattering phase shifts evaluated at the on-shell relative momenta derived 
from the spectrum. In EigTI^we give illustrative examples of the quality of the results for the spin-triplet and spin- 
singlet channels in various irreps that couple predominantly to a given partial wave. The bands are fits to the effective 
range expansion (ERE) to different orders, i.e. 

cot (5^ = -b -b yPeq'^ + 0{q^), (5) 

ae 2 4! 

where ai, r^, and Pi are the scattering length, effective range, and shape parameter for i = 0 and the corresponding 
parameters of the ERE for i > 0. Results for LO {q^) fits are denoted by yellow bands, fits to NLO are blue, and 
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channel 

1 / [a,ml^+^) 

n 


ERE order 

channel 


ERE order 


0.088(10) 

4.47(31) 

- 


"Ho 

0.234(75) 

q^ 


0.056(24) 

5.45(76) 

-202(114) 



0.237(92) 

q° 

"S'! 

0.094(06) 

4.27(21) 

- 

q^ 


-0.317(49) 

q° 


0.071(09) 

4.51(24) 

-139(38) 

q^ 


0.146(22) 

q° 

"U 3 

-0.046(20) 

- 

- 


^D2 

-0.047(12) 

q° 


-0.082(51) 

-0.27(34) 

- 

q^ 





TABLE I: Effective range parameters for the various scattering channels determined from fits to the orders indicated in the 
last column. 


NNLO (q^) fits are red. The ERE parameters determined from these fits are listed in Table The dashed vertical 
line represents the t-channel cut, above which the ERE is expected to break down. While the ERE becomes formally 
unjustifiable past this cut, the Liischer formalism holds for all energies below the NNtt threshold, well above the 
energies considered. 

We obtain results significantly different from zero for various channels, including P and D channels. One beautiful 
illustration of the power of the operators and the finite volume formalism being used is that of ^P 2 in Fig. 0 (top- 
middle and top-right). In this channel the spectrum was determined with two irreps that have overlap with the same 
partial wave using two different volumes. The consistency of the extracted phase shifts clearly demonstrates that the 
generalization of the Liischer formalism for NN-systems is working [53] . We hnd that over the kinematic range of our 
calculations, cot Sap^ is consistent with a constant, even above the t-channel cut. Furthemore, we find no evidence 
of the t-channel cut playing an important role for any of the channels studied at these values of the quark masses. 

Currently, we ignore partial wave mixing, an issue that will be addressed in subsequent publications. This formidable 
challenge has only been addressed in two-meson calculations [Ml IM IlSHIHj- However, some evidence that the 
mixing from higher partial waves is small in at least one channel can be obtained by investigating the results for the 
^P 2 channel in Fig. Again, for this channel we have two cubic irreps, and E~, for which the lowest partial wave 
is ^P 2 - Any differences between the two irreps must arise from mixing with higher partial waves. We find that the 
two cubic irreps give completely consistent results, indicating that this mixing must be too small to resolve within 
our error bars. Additionally, the first contribution from mixing in the cubic irrep, having lowest partial wave ^Pi, 
comes from the partial wave. We obtain information on the strength of the phase shift independently using 
the A 2 cubic irrep, and find it to be extremely small, thus it likely does not contribute to mixing in the Tf cubic 
irrep. 


A. Bound states from the Effective Range Expansion 

Assuming the results are within the range of convergence of the ERE, the infinite volume bound state energies may 
be determined by solving for poles in the derived scattering amplitude. For a single channel scattering amplitude, 
these satisfy 


qp cot 6{qB) = iqs , (6) 

which has a solution below threshold for qp = inp- In both S'-wave channels, the ERE expansion appears to be 
converging well (top-left and middle-left figures in Fig. with small 0{q^) corrections where the results exist below 
the t-channel cut. Using this method and the spectrum obtained in this work, we find deeply bound states in both 
the ^S'o and channels, with binding energies of 

Bas, = 21.8(^|;?)(t°i) MeV, Bap, = 30.7(t2-4)(+o.5) MeV. (7) 

The first uncertainties are our fitting statistical and systematic uncertainties combined in quadrature and the second 
is an estimate of systematic uncertainties arising from the truncation of the ERE. In addition, in the channel we 
find a second pole near threshold, corresponding to 

=3.3(J:°)(°:®)MeV. (8) 

Corresponding to each of these poles we find finite volume states whose energies are consistent with the expected 
exponential volume dependence associated with bound states. However, with only two volumes we cannot definitively 
state whether the volume dependence is exponential or polynomial. With this precision, it is unclear whether this 
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L 

unitarity bound 

T+{^Si) 





-3.786/(ML2) 

AEo 

AEi 

AEo 

AEi 

24 

-7.41 

-30.4(2.4)(5.1) 

-21.4(1.0)(0.5) 

-20.2(2.1)(1.5) 

- 

32 

-4.17 

-28.1(1.8)(2.4) 

-9.95(.99)b42) 

-17.3(1.7)(2.3) 

-8.35(.99)(.48) 


TABLE II: Negative shifted energy states in ^Si and channels compared with the unitarity bound. These energy levels 
have been converted to MeV using the lattice spacing b — 0.145 fm. No scale setting uncertainty is assigned as we are just 
comparing the energy levels to the unitarity bound. 


shallow bound state corresponds to a true bound state or a near-threshold scattering state. Improved analysis 
techniques, such as a full basis of interpolating helds in momentum space, as has been successfully used in the two- 
meson systems [551 UHl H5M5] . or additional statistics may be necessary to fully settle this matter. We do note, 
however, that all of the negative energy finite-volume shifts are larger by 4 or more standard deviations than the 
corresponding finite-volume energy shift for a bound state at threshold, AE ^ |711|7S], as shown in Table [n| 

This indicates that the interactions producing these states are more attractive than those at unitarity, and the states 
should therefore correspond to true bound states in the infinite volume limit. Further evidence of multiple bound 
states are presented in the next section. 


B. Evidence for multiple negatively shifted energy states 

We further elaborate on the plausibility of finding two negatively shifted energy states in both S'-wave channels. 
The large negatively shifted energy levels, determined with the local operators, are consistent with those in Ref. [57], 
which also used local operators.^ The state closer to threshold (and additionally, the negative energy state near 
threshold in the channel) has strong overlap onto the non-local NN interpolating field, and has not been found 
in previous works. In Fig. [^ we plot the effective mass of the ratio correlation functions of the non-local two-nucleon 
interpolating fields divided by the local two-nucleon interpolating fields 

RyTvW = C'h5^oW/C'm=o(0, C'hW = 

n 

Both two-nucleon interpolating fields couple to the same tower of states, = 2mN -I- A)) and only differ in the 
relative size of their overlap factors, A))(|r|). In the spin triplet, T(*" channel, there is a clear plateau inconsistent with 
0 for a relatively large time in fm (recall that each time step on these ensembles correspond to a = 0.1453(16) fm [57] b 
indicating that the two correlation functions each give statistically distinct plateaus. For the spin singlet, channel, 
a clear plateau is only observed on the larger volume. 

In the long-time limit, the effective masses for all interpolating operators in this channel must asymptote to the 
ground state of the system. However, if the non-local interpolating field couples strongly to the state with small 
negative energy shift, and weakly to the state with large negative energy shift, and conversely for the local interpolating 
field, then at intermediate times the dominant contribution to each interpolating field will come from the two different 
states respectively. This would manifest in the non-local interpolating held having an effective mass plateau at 
the smaller negative energy shift (the hrst excited state) at intermediate times. With enough statistics, one would 
eventually observe the effective mass “collapse” to the ground state and plateau again at the larger negative energy 
shift (the ground state). Note that this is also true of all excited two-nucleon elastic scattering states. Because relative 
momentum is not a good quantum number, the non-zero relative momentum projections serve only to enhance the 
overlap of the interpolating operators with excited states relative to the ground state. We hnd this to work quite 
well due to the existence of distinct plateaus at intermediate times for the different relative momentum operators, 
however, in the very long time limit we expect all of these correlation functions to collapse onto a single ground state 
plateau. 

This explanation for the two states with negative energies is consistent with the intuition that the wavefunction 
for a shallow bound state is more extended in space, and thus will have poor overlap with an interpolating field 
involving only two nucleons at the same point (a smeared out delta function in our case). Conversely, the more deeply 


^ Note, our results for the binding energies do differ somewhat from those in Ref. El due in part to different strategies to determine the 
infinite volume binding energies from the finite volume energy levels. 
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FIG. 5: Differences of effective masses using displaced vs. local sources for the lowest energy levels in the spin singlet, Af and 
spin triplet,channels on each volume. In all cases except for the spin singlet, Af channel with L/h = 24, we find signihcantly 
distinct energy levels using the two different types of sources. 


bound state would have a much more compact wave function and would thus overlap more poorly with a non-local 
interpolating field. We finally emphasize that these findings are consistent with the two poles found using the ERE 
of our phase shift, even if we disregard the finite volume negative energy results from either the local or the non-local 
operators. We therefore find the most plausible explanation of these results to be the existence of two distinct bound 
states in the channel. 


IV. SUMMARY 

This work presents the implementation of new two-nucleon interpolating fields which allow, for the first time, a 
robust determination of £ > 0 scattering phase shifts in the NN sector, which we have calculated in this exploratory 
work using unphysically heavy quark masses. Further, this improved basis of interpolating operators are sensitive to 
additional states in the S-wave spectrum that were not found using only local operators and greater statistics. This 
has been made possible by three previously unexploited tools. First was the development of displaced two-nucleon 
interpolating sources. These are necessary to have appreciable overlap with partial waves beyond the S-wave as the 
£ 0 orbital wavefunctions are zero at the origin. Second was the use of momentum space sink operators that were 

not restricted to the simplest cubic irreps. Finally, we applied the formalism for two-nucleon systems in a finite 
volume [69) . ignoring higher partial wave mixing, with notable success for the ^P 2 channel. This work represents the 
first crucial step towards the study of more challenging systems such as three-neutron interactions and the S ^ P 
wave parity violating pp interaction. 
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